#################
# Lebovic and Voeten Replication
# with Causal Endogeneity Model
#################

rm(list=ls())
set.seed(123456)
library(foreign)
library(lme4)
library(mvtnorm)
dat<-read.dta("jprworkdatanewMOD.dta")

# the Lebovic and Voeten replication model
col1<-lmer(BIPOP~l_BIPOP+PUBRES+d_HRIGHTS+l_HRIGHTS+d_CIVIL+l_CIVIL+l_GDPPOP+l_LNPOP+l_USAGREE+WAR+CAPAB+linear+quad+(1 | CCODE), data=dat, REML=F)
options(scipen=100)
summary(col1)

# the state dependence model
causend<-lmer(BIPOP~l_BIPOP+PUBRES+RESxBI+d_HRIGHTS+l_HRIGHTS+d_CIVIL+l_CIVIL+l_GDPPOP+l_LNPOP+l_USAGREE+WAR+CAPAB+linear+quad+(1 | CCODE), data=dat, REML=F)
options(scipen=100)
summary(causend)

# prepare for marginal effects calculations
beta.draws<-rmvnorm(1000, mean=causend@beta, sigma=as.matrix(vcov(causend)))

# instantaneous effects of a CHR Resolution
y.seq<-seq(from=0, to=7.25, by=0.1)
me<-matrix(data=NA, ncol=3, nrow=length(y.seq))

for(i in 1:length(y.seq)){
  
  y.test<-y.seq[i]
  me.draws<-beta.draws[,3]+beta.draws[,4]*y.test
  me[i,]<-quantile(me.draws, probs=c(0.025, 0.5, 0.975))
  
  
}

# plot instantaneous marginal effects
pdf("instant_mfx_applied_RE_r.pdf")
par(mar=c(5,5,4,2))
plot(NA, ylim=c(-0.6, 1), xlim=c(0,7.5), xaxt="n", main=c("Instantaneous Marginal Effect","of UNCHR Resolution on Aggregate Bilateral Aid"), xlab=expression(ln(Aid + 1)[it-1]), ylab=expression(partialdiff*ln(Aid + 1)[it]/partialdiff*Resolution[i(t-1)]))
axis(side=1, at=0:7)
hist.back<-hist(dat$l_BIPOP, plot=FALSE)
lines(hist.back, border=hsv(0, 0, 0.5, alpha=1), col=hsv(0, 0, 0.5, alpha=0.25), freq=FALSE, alpha=1)

lines(me[,2]~y.seq)
lines(me[,1]~y.seq, lty=2)
lines(me[,3]~y.seq, lty=2)
abline(h=0, lty=3)

legend("topleft", legend=c("marginal effect", "95% confidence interval"), lty=c(1,2))
text(x=5.75, y=-0.6, "histogram indicates data density", cex=0.75)
dev.off()




# long term marginal effects of a CHR resolution
z<-seq(from=0, to=3.25, by=0.05)
ss.orig<-c()
ss.hi<-c()
ss.lo<-c()
ss.med<-c()
for(i in 1:length(z)){
  
  x.lo<-0
  x.hi<-1
  
  ss.orig[i]<-mean((z[i]+beta.draws[,3]*x.lo)/(1-beta.draws[,2]-beta.draws[,4]*x.lo))
  ss.draws<-((z[i]+beta.draws[,3]*x.hi)/(1-beta.draws[,2]-beta.draws[,4]*x.hi))-((z[i]+beta.draws[,3]*x.lo)/(1-beta.draws[,2]-beta.draws[,4]*x.lo))
  ss.quant<-quantile(ss.draws, probs=c(0.025, 0.5, 0.975))
  ss.hi[i]<-ss.quant[3]
  ss.lo[i]<-ss.quant[1]
  ss.med[i]<-ss.quant[2]
  
}

summary(ss.orig)

pdf("longtermME_applied_r.eps")
plot(ss.med~ss.orig, type="l", ylim=c(-3, 4), main=c("Long Term Marginal Effect", "of UNCHR Resolution on Bilateral Aid"), xlab="ln(Bilateral Aid + 1) steady state, no resolution", ylab="change in ln(Bilateral Aid + 1) steady state after resolution")
lines(ss.hi~ss.orig, lty=2)
lines(ss.lo~ss.orig, lty=2)
abline(h=0, lty=3)
legend("bottomright", lty=c(1,2), legend=c("marginal effect", "95% confidence interval"))
dev.off()






ss.hi<-c()
ss.lo<-c()
ss.med<-c()
ss.draws<-matrix(data=NA, ncol=2, nrow=1000)


# long term effect of a human rights violation
x<-c(0,1)
for(i in 1:2){
  ss.draws[,i]<-(beta.draws[,6]*4)/(1-beta.draws[,2]-beta.draws[,4]*x[i])
  ss.quant<-quantile(ss.draws[,i], probs=c(0.025, 0.5, 0.975))
  ss.hi[i]<-ss.quant[3]
  ss.lo[i]<-ss.quant[1]
  ss.med[i]<-ss.quant[2]
}

library(gplots)
pdf("longtermME_Z_applied_r.eps")
par(mar=c(4, 5, 4, 2))
plotCI(x=x, y=ss.med, ui=ss.hi, li=ss.lo, ylim=c(-0.85, 0.25), xlim=c(-0.5, 1.5), xaxt="n", yaxt="n", main=c("Long Term Marginal Effect", "of Physical Integrity Abuse on Bilateral Aid"), ylab="", xlab="")
axis(1, at=0:1, lab=c("w/o UNCHR Resolution","w/ UNCHR Resolution"))
axis(2, at=seq(from=-0.8, to =0.2, by=0.2))
mtext(text="change in ln(Bilateral Aid + 1) steady state,", side=2, line=4)
mtext(text="abuse from min to max", side=2, line=3)
abline(h=0, lty=3)
dev.off()









# plot long term trajectories for the modified LV model

x<-matrix(dat=c(1,1,0,0,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
y<-rbind(1,x%*%t(beta.draws)-0.75)
i<-2
k<-2

while(round(y[i,1], 7)!=round(y[i-1,1],7)){
  y<-rbind(y,NA)
  i<-i+1
  k<-k+1
  for(j in 1:1000){
    x<-matrix(dat=c(1,y[i-1,j],0,y[i-1,j]*0,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
    y[i,j]<-x%*%t(t(beta.draws[j,]))-0.75
  }
  
}


y<-rbind(y,NA)
i<-i+1
for(j in 1:1000){
  x<-matrix(dat=c(1,y[i-1,j],1,y[i-1,j]*1,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
  y[i,j]<-x%*%t(t(beta.draws[j,]))-0.75
}

pen<-i+9

while(i<pen){
  i<-i+1
  y<-rbind(y,NA)
  for(j in 1:1000){
    x<-matrix(dat=c(1,y[i-1,j],1,y[i-1,j]*1,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
    y[i,j]<-x%*%t(t(beta.draws[j,]))-0.75
  }
}

while(round(y[i,1], 7)!=round(y[i-1,1],7)){
  y<-rbind(y,NA)
  i<-i+1
  for(j in 1:1000){
    x<-matrix(dat=c(1,y[i-1,j],0,y[i-1,j]*0,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
    y[i,j]<-x%*%t(t(beta.draws[j,]))-0.75
  }
  
}

end.point<-(k-10)
t<-0:(i-end.point)
plot(y[end.point:i]~t, type="l", xlab="time")
abline(v=10, lty=2)
abline(v=20, lty=2)

traj.quant<-apply(X=y, MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))

plot(traj.quant[2,end.point:i]~t, type="l", xlab="time")
lines(traj.quant[1,end.point:i]~t, lty=2)
lines(traj.quant[3,end.point:i]~t, lty=2)

library(gplots)
pdf("smallaidtrajectory.eps")
plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=traj.quant[1,end.point:i], ylab="ln(aggregate bilateral aid PC + 1)", xlab="time", ylim=c(.41, 1.81), main = c("simulated aid trajectory with 90% CI error bars"))
rect(11,-1,20,4,density=NA, col="grey85")
plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=traj.quant[1,end.point:i], add=T)
box()
abline(v=11, lty=2, lwd=2)
abline(v=20, lty=2, lwd=2)
legend("bottomright", legend=c(" median aid prediction", " condemnation period"), pt.cex=c(1,2), pch=c(1,22), pt.bg="gray85")
dev.off()

traj.quant[2, (end.point+10):(end.point+20)]
# percentage change in aid
1 - exp(traj.quant[2, (end.point+20)]) / exp(traj.quant[2, (end.point+10)])






x<-matrix(dat=c(1,1,0,0,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
y<-rbind(1,x%*%t(beta.draws)+0.5)
i<-2
k<-2

while(round(y[i,1], 7)!=round(y[i-1,1],7)){
  y<-rbind(y,NA)
  i<-i+1
  k<-k+1
  for(j in 1:1000){
    x<-matrix(dat=c(1,y[i-1,j],0,y[i-1,j]*0,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
    y[i,j]<-x%*%t(t(beta.draws[j,]))+0.5
  }
  
}


y<-rbind(y,NA)
i<-i+1
for(j in 1:1000){
  x<-matrix(dat=c(1,y[i-1,j],1,y[i-1,j]*1,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
  y[i,j]<-x%*%t(t(beta.draws[j,]))+0.5
}

pen<-i+9

while(i<pen){
  i<-i+1
  y<-rbind(y,NA)
  for(j in 1:1000){
    x<-matrix(dat=c(1,y[i-1,j],1,y[i-1,j]*1,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
    y[i,j]<-x%*%t(t(beta.draws[j,]))+0.5
  }
}

while(round(y[i,1], 7)!=round(y[i-1,1],7)){
  y<-rbind(y,NA)
  i<-i+1
  for(j in 1:1000){
    x<-matrix(dat=c(1,y[i-1,j],0,y[i-1,j]*0,0,2.785247,0,4.692341,6.728779,15.8313, .3541963,0,0.0233268,0,0), nrow=1)
    y[i,j]<-x%*%t(t(beta.draws[j,]))+0.5
  }
  
}

end.point<-(k-10)

t<-0:(i-end.point)
plot(y[end.point:i,1]~t, type="l", xlab="time")
abline(v=10, lty=2)
abline(v=20, lty=2)

traj.quant<-apply(X=y, MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))


library(gplots)
pdf("largeaidtrajectory.eps")
plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=traj.quant[1,end.point:i], ylab="ln(aggregate bilateral aid PC + 1)", xlab="time", main = c("simulated aid trajectory with 90% CI error bars"), ylim=c(3.99, 5.41))
rect(11,0,20,6,density=NA, col="grey85")
plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=traj.quant[1,end.point:i], add=T)
box()
abline(v=11, lty=2, lwd=2)
abline(v=20, lty=2, lwd=2)
legend("topright", legend=c(" median aid prediction", " condemnation period"), pt.cex=c(1,2), pch=c(1,22), pt.bg="gray85")
dev.off()

traj.quant[2, (end.point+10):(end.point+20)]
# percentage change in aid
exp(traj.quant[2, (end.point+20)]) / exp(traj.quant[2, (end.point+10)]) - 1

















### plot instantaneous MEs, Neilsen Model 3


rm(list=ls())
library(MASS)
require(foreign)

dat.neilsen<-read.dta(file="neilsen_3_R_figure_data.dta")

neilsen.beta<-read.delim(file="xttobit_beta.txt", header=T)
neilsen.beta<-as.matrix(neilsen.beta[2:31])
names(neilsen.beta)<-NULL

neilsen.vcv<-read.delim(file="xttobit_VCV.txt")
neilsen.vcv<-as.matrix(neilsen.vcv[1:30,2:31])

beta.draws<-mvrnorm(n=1000, mu=neilsen.beta, Sigma=neilsen.vcv)


# instantaneous effects of a CHR Resolution
par(mar=c(5,5,4,2))
y.seq<-seq(from=0, to=13, by=0.1)
me<-matrix(data=NA, ncol=3, nrow=length(y.seq))

for(i in 1:length(y.seq)){
  
  y.test<-y.seq[i]
  me.draws<-beta.draws[,1]+beta.draws[,2]*y.test
  me[i,]<-quantile(me.draws, probs=c(0.025, 0.5, 0.975))
  
  
}

pdf("neilsen_instant.pdf")
# plot instantaneous marginal effects
plot(NA, ylim=c(-2, 2.5), xlim=c(0,13), xaxt="n", main=c("Instantaneous Marginal Effect","of UNCHR Resolution on Bilateral Economic Aid"), xlab=expression(ln("Bilateral Economic Aid + 1")[it-1]), ylab=expression(partialdiff*ln("Bilateral Economic Aid + 1")[it]/partialdiff*Resolution[i(t-1)]))
axis(side=1, at=seq(from=0, to=12, by=2))
hist.back<-hist(dat.neilsen$l_lneconaidpc, plot=FALSE)
hist.back$density <- hist.back$density*3
lines(hist.back, border=hsv(0, 0, 0.5, alpha=1), col=hsv(0, 0, 0.5, alpha=0.25), freq=FALSE, alpha=1)

lines(me[,2]~y.seq)
lines(me[,1]~y.seq, lty=2)
lines(me[,3]~y.seq, lty=2)
abline(h=0, lty=3)

legend("topleft", legend=c("marginal effect", "95% confidence interval"), lty=c(1,2))
text(x=10, y=-2, "histogram indicates 3*(data density)", cex=0.75)
dev.off()


### plot long term trajectories, Neilsen model 3

x<-matrix(dat=c(0, 2*0, 3.880258, 0, 0*0, 0, 0*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, 2, 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
y<-rbind(1,x%*%t(beta.draws)+8.8)
i<-2
k<-2

while(round(y[i,1], 7)!=round(y[i-1,1],7)){
  y<-rbind(y,NA)
  i<-i+1
  k<-k+1
  for(j in 1:1000){

    x<-matrix(dat=c(0, y[i-1,j]*0, 3.880258, 0, 0*0, 0, 0*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
    y[i,j]<-x%*%t(t(beta.draws[j,]))+8.8
    #print(y[i,1], quote=F)
    
  }
  
}


y<-rbind(y,NA)
i<-i+1
for(j in 1:1000){
  
  x<-matrix(dat=c(1, y[i-1,j]*1, 3.880258, 0, 0*1, 0, 0*1, .3926666, .3926666*1, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
  y[i,j]<-x%*%t(t(beta.draws[j,]))+8.8
  #print(y[i,1], quote=F)
  
}

pen<-i+9

while(i<pen){
  
  y<-rbind(y,NA)
  i<-i+1
  for(j in 1:1000){
    
    x<-matrix(dat=c(1, y[i-1,j]*1, 3.880258, 0, 0*1, 0, 0*1, .3926666, .3926666*1, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
    
    y[i,j]<-x%*%t(t(beta.draws[j,]))+8.8
    #print(y[i,1], quote=F)
    
  }
  
}


while(round(y[i,1], 7)!=round(y[i-1,1],7)){
  y<-rbind(y,NA)
  i<-i+1
  for(j in 1:1000){
    
    x<-matrix(dat=c(0, y[i-1,j]*0, 3.880258, 0, 0*0, 0, 0*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
    y[i,j]<-x%*%t(t(beta.draws[j,]))+8.8
    #print(y[i,1], quote=F)
    
  }
  
}

end.point<-(k-10)
t<-0:(i-end.point)
plot(y[end.point:i,1]~t, type="l", xlab="time")
abline(v=10, lty=2)
abline(v=20, lty=2)


traj.quant<-apply(X=y, MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))


pdf("dyadaidtrajectory.eps")
library(gplots)
plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=traj.quant[1,end.point:i], ylab="ln(bilateral aid PC + 1)", xlab="time", main = c("simulated aid trajectory with 90% CI error bars"), ylim=c(5.78, 10.02))
rect(11,0,20,15,density=NA, col="grey85")
plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=traj.quant[1,end.point:i], add=T)
box()
abline(v=11, lty=2, lwd=2)
abline(v=20, lty=2, lwd=2)
legend("topright", legend=c(" median aid prediction", " condemnation period"), pt.cex=c(1,2), pch=c(1,22), pt.bg="gray85")
dev.off()


traj.quant[2, (end.point+10):(end.point+20)]
# percentage change in aid
exp(traj.quant[2, (end.point+20)]) / exp(traj.quant[2, (end.point+10)]) - 1





x<-matrix(dat=c(0, 2*0, 3.880258, 0, 0*0, 0, 0*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, 2, 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
y<-rbind(1,x%*%t(beta.draws)+6)
i<-2
k<-2

while(round(y[i,1], 7)!=round(y[i-1,1],7)){
  y<-rbind(y,NA)
  i<-i+1
  k<-k+1
  for(j in 1:1000){
    
    x<-matrix(dat=c(0, y[i-1,j]*0, 3.880258, 0, 0*0, 0, 0*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
    y[i,j]<-x%*%t(t(beta.draws[j,]))+6
    #print(y[i,1], quote=F)
    
  }
  
}


y<-rbind(y,NA)
i<-i+1
for(j in 1:1000){
  
  x<-matrix(dat=c(1, y[i-1,j]*1, 3.880258, 0, 0*1, 0, 0*1, .3926666, .3926666*1, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
  y[i,j]<-x%*%t(t(beta.draws[j,]))+6
  #print(y[i,1], quote=F)
  
}

pen<-i+9

while(i<pen){
  
  y<-rbind(y,NA)
  i<-i+1
  for(j in 1:1000){
    
    x<-matrix(dat=c(1, y[i-1,j]*1, 3.880258, 0, 0*1, 0, 0*1, .3926666, .3926666*1, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
    y[i,j]<-x%*%t(t(beta.draws[j,]))+6
    #print(y[i,1], quote=F)
    
  }
  
}


while(round(y[i,1], 7)!=round(y[i-1,1],7)){
  y<-rbind(y,NA)
  i<-i+1
  for(j in 1:1000){
    
    x<-matrix(dat=c(0, y[i-1,j]*0, 3.880258, 0, 0*0, 0, 0*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
    y[i,j]<-x%*%t(t(beta.draws[j,]))+6
    #print(y[i,1], quote=F)
    
  }
  
}

end.point<-(k-10)
t<-0:(i-end.point)
plot(y[end.point:i,1]~t, type="l", xlab="time")
abline(v=10, lty=2)
abline(v=20, lty=2)


traj.quant<-apply(X=y, MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))


library(gplots)
pdf("dyadaidtrajectory_small.eps")
plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=pmax(traj.quant[1,end.point:i],0), ylab="ln(bilateral aid PC + 1)", xlab="time", main = c("simulated aid trajectory with 90% CI error bars"))
rect(11,-1,20,15,density=NA, col="grey85")
plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=pmax(traj.quant[1,end.point:i],0), add=T)
box()
abline(v=11, lty=2, lwd=2)
abline(v=20, lty=2, lwd=2)
legend("bottomright", legend=c(" median aid prediction", " condemnation period"), pt.cex=c(1,2), pch=c(1,22), pt.bg="gray85")
dev.off()


traj.quant[2, (end.point+10):(end.point+20)]
# percentage change in aid
1 - exp(traj.quant[2, (end.point+20)]) / exp(traj.quant[2, (end.point+10)])








x<-matrix(dat=c(0, 2*0, 3.880258, 0, 0*0, 1, 1*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, 2, 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
y<-rbind(1,x%*%t(beta.draws)+6)
i<-2
k<-2

while(round(y[i,1], 7)!=round(y[i-1,1],7)){
  y<-rbind(y,NA)
  i<-i+1
  k<-k+1
  for(j in 1:1000){
    
    x<-matrix(dat=c(0, y[i-1,j]*0, 3.880258, 0, 0*0, 1, 1*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
    y[i,j]<-x%*%t(t(beta.draws[j,]))+6
    #print(y[i,1], quote=F)
    
  }
  
}


y<-rbind(y,NA)
i<-i+1
for(j in 1:1000){
  
  x<-matrix(dat=c(1, y[i-1,j]*1, 3.880258, 0, 0*1, 1, 1*1, .3926666, .3926666*1, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
  y[i,j]<-x%*%t(t(beta.draws[j,]))+6
  #print(y[i,1], quote=F)
  
}

pen<-i+9

while(i<pen){
  
  y<-rbind(y,NA)
  i<-i+1
  for(j in 1:1000){
    
    x<-matrix(dat=c(1, y[i-1,j]*1, 3.880258, 0, 0*1, 1, 1*1, .3926666, .3926666*1, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
    y[i,j]<-x%*%t(t(beta.draws[j,]))+6
    #print(y[i,1], quote=F)
    
  }
  
}


while(round(y[i,1], 7)!=round(y[i-1,1],7)){
  y<-rbind(y,NA)
  i<-i+1
  for(j in 1:1000){
    
    x<-matrix(dat=c(0, y[i-1,j]*0, 3.880258, 0, 0*0, 1, 1*0, .3926666, .3926666*0, 4.00579,  1.257437, .7319112,  7.442243, -.5143044, y[i-1,j], 15.93863, 7.850649, 9.193985, 15.96552, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), nrow=1)
    y[i,j]<-x%*%t(t(beta.draws[j,]))+6
    #print(y[i,1], quote=F)
    
  }
  
}

end.point<-(k-10)
t<-0:(i-end.point)
plot(y[end.point:i,1]~t, type="l", xlab="time")
abline(v=10, lty=2)
abline(v=20, lty=2)


traj.quant<-apply(X=y, MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))


pdf("dyadaidtrajectory_alliance.eps")
library(gplots)
plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=pmax(traj.quant[1,end.point:i],0), ylab="ln(bilateral aid PC + 1)", xlab="time", main = c("simulated aid trajectory with 90% CI error bars"))
rect(11,-1,20,15,density=NA, col="grey85")
plotCI(x=t, y=traj.quant[2,end.point:i], ui=traj.quant[3,end.point:i], li=pmax(traj.quant[1,end.point:i],0), add=T)
box()
abline(v=11, lty=2, lwd=2)
abline(v=20, lty=2, lwd=2)
legend("topright", legend=c(" median aid prediction", " condemnation period"), pt.cex=c(1,2), pch=c(1,22), pt.bg="gray85")
dev.off()


traj.quant[2, (end.point+10):(end.point+20)]
# percentage change in aid
1 - exp(traj.quant[2, (end.point+20)]) / exp(traj.quant[2, (end.point+10)])










### plot instantaneous MEs, Murdie/Davis Model 2


rm(list=ls())
library(MASS)
require(foreign)

dat.neilsen <- read.dta(file="ngo-shaming.dta")

neilsen.beta<-read.delim(file="tobit_beta_ngo.txt", header=T)
neilsen.beta<-as.matrix(neilsen.beta[2:26])
names(neilsen.beta)<-NULL

neilsen.vcv<-read.delim(file="tobit_VCV_ngo.txt")
neilsen.vcv<-as.matrix(neilsen.vcv[1:25,2:26])

beta.draws<-mvrnorm(n=1000, mu=neilsen.beta, Sigma=neilsen.vcv)


# instantaneous effects of NGO Shaming
par(mar=c(5,5,4,2))
y.seq<-seq(from=0, to=13, by=0.1)
me<-matrix(data=NA, ncol=3, nrow=length(y.seq))

for(i in 1:length(y.seq)){
  
  y.test<-y.seq[i]
  me.draws<-beta.draws[,1]+beta.draws[,2]*y.test
  me[i,]<-quantile(me.draws, probs=c(0.025, 0.5, 0.975))
  
  
}

pdf("ngo_instant.pdf")
# plot instantaneous marginal effects
plot(NA, ylim=c(-0.25, 0.25), xlim=c(0,13), xaxt="n", main=c("Instantaneous Marginal Effect","of NGO Shaming on Bilateral Economic Aid"), xlab=expression(ln("Bilateral Economic Aid + 1")[it-1]), ylab=expression(partialdiff*ln("Bilateral Economic Aid + 1")[it]/partialdiff*Shaming[i(t-1)]))
axis(side=1, at=seq(from=0, to=12, by=2))
hist.back<-hist(dat.neilsen$l_lneconaidpc, plot=FALSE)
hist.back$density <- hist.back$density*0.3
lines(hist.back, border=hsv(0, 0, 0.5, alpha=1), col=hsv(0, 0, 0.5, alpha=0.25), freq=FALSE, alpha=1)

lines(me[,2]~y.seq)
lines(me[,1]~y.seq, lty=2)
lines(me[,3]~y.seq, lty=2)
abline(h=0, lty=3)

legend("topright", legend=c("marginal effect", "95% confidence interval"), lty=c(1,2))
text(x=10, y=-2, "histogram indicates 3*(data density)", cex=0.75)
dev.off()













rm(list=ls())
library(MASS)
require(foreign)

dat.neilsen <- read.dta(file="ngo-shaming.dta")

neilsen.beta<-read.delim(file="tobit_beta_ngo2.txt", header=T)
neilsen.beta<-as.matrix(neilsen.beta[2:26])
names(neilsen.beta)<-NULL

neilsen.vcv<-read.delim(file="tobit_VCV_ngo2.txt")
neilsen.vcv<-as.matrix(neilsen.vcv[1:25,2:26])

beta.draws<-mvrnorm(n=1000, mu=neilsen.beta, Sigma=neilsen.vcv)


# instantaneous effects of NGO Shaming
par(mar=c(5,5,4,2))
y.seq<-seq(from=0, to=13, by=0.1)
me<-matrix(data=NA, ncol=3, nrow=length(y.seq))

for(i in 1:length(y.seq)){
  
  y.test<-y.seq[i]
  me.draws<-beta.draws[,1]+beta.draws[,2]*y.test
  me[i,]<-quantile(me.draws, probs=c(0.025, 0.5, 0.975))
  
  
}

pdf(file="ngo_dum_instant.pdf")
# plot instantaneous marginal effects
plot(NA, ylim=c(-1.5, 1.5), xlim=c(0,13), xaxt="n", main=c(""), xlab=expression(ln("Bilateral Economic Aid + 1")[it-1]), ylab=expression(partialdiff*ln("Bilateral Economic Aid + 1")[it]/partialdiff*Shaming[i(t-1)]))
mtext(line=0.75, expression(paste("NGO Shaming">=1~"on Bilateral Economic Aid")), cex=1.1)
mtext(line=2, "Instantaneous Marginal Effect of", cex=1.1)
axis(side=1, at=seq(from=0, to=12, by=2))
hist.back<-hist(dat.neilsen$l_lneconaidpc, plot=FALSE)
hist.back$density <- hist.back$density*2
lines(hist.back, border=hsv(0, 0, 0.5, alpha=1), col=hsv(0, 0, 0.5, alpha=0.25), freq=FALSE, alpha=1)

lines(me[,2]~y.seq)
lines(me[,1]~y.seq, lty=2)
lines(me[,3]~y.seq, lty=2)
abline(h=0, lty=3)

legend("topright", legend=c("marginal effect", "95% confidence interval"), lty=c(1,2))
text(x=10, y=-2, "histogram indicates 3*(data density)", cex=0.75)

dev.off()

